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RESEARCH LABORATORIES FOR THE ENGINEERING SCIENCES 


Members of the faculty who teach at the undergraduate and graduate levels and a number of 
professional engineers and scientists whose primary activity is research generate and conduct the 
investigations that make up the school’s research program. The School of Engineering and Applied Science 
of the University of Virginia believes that research goes hand in hand with teaching. Early in the 
development of its graduate training program, the School recognized that men and women engaged in 
research should be as free as possible of the administrative duties involved in sponsored research. In 1959, 
therefore, the Research Laboratories for the Engineering Sciences (RLES) was established and assigned the 
administrative responsibility for such research within the School. 


The director of RLES—himself a faculty member and researcher—maintains familiarity with the 
support requirements of the research under way. He is aided by an Academic Advisory Committee made up 
of a faculty representative from each academic department of the School. This Committee serves to inform 
RLES of the needs and perspectives of the research program. 


In addition to administrative support, RLES is charged with providing certain technical assistance. 
Because it is not practical for each department to become self-sufficient in all phases of the supporting 
technology essential to present-day research, RLES makes services available through the following support 
groups: Machine Shop, Instrumentation, Facilities Services, Publications (including photographic facilities), 
and Computer Terminal Maintenance. 
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1.0 SUMMARY 


This research and development effort has obtained the following 
results: 

* Software suitable for the construction of finite element gravity 
disturbance fields has been developed. 

* Software suitable for calculation of the gravity disturbance model 
within finite element fields has been developed. 

¢ Example finite element gravity disturbance fields have been 


developed and tested. 
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2.0 PREFACE 


This report constitutes the final report of Phase II of Contract 
No. DAAG-53-76-C-0067 performed by the University of Virginia for the U. S. 
Army Engineer Topographic Laboratories, Fort Belvoir, Virginia, under the 
sponsorship of the Defense Mapping Agency. The Phase II effort is con- 
cerned with representing local fine structure gravity by locally valid 
(finite element) gravity functions. The Phase I effort, documented in 
Reference (1), deals with replacing globally valid gravity representations 
(such as spherical harmonic series) by an equivalent family of locally 
valid gravity functions. The Phase I effort yielded a macro gravity model 
suitable for high speed satellite orbit integration. The Phase II effort 
(documented herein) yielded fine structure gravity representations for 
efficient disturbance acceleration calculation (to use, for example, 


in inertial guidance equations). 


The authors appreciate the technical inputs of Mr. L. A. Gambino 
(technical monitor), Dr. R. W. Ballew, Mr. H. W. Howard, and Mr. B. L. Decker. 
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3.0 INTRODUCTION 


The finite element concept may be considered as an extension to 
three dimensions of one dimensional piecewise approximation techniques. 
Some quantity that is functionally related in a "complicated" fashion to 
position over a “large” interval may be replaced for much "smaller" 
intervals of position with much "less complicated" approximating functions. 
Thus the complicated function, which will normally be “more expensive" to 
evaluate, will have been replaced by a set of locally valid functions 


which will normally be "less expensive" to evaluate. 


The idea may easily be extended to three dimensions. A quantity 
which is functionally related to position over some large region may be 
replaced by a set of approximations which are valid only in small volumes 
of the large region, but which are easier (faster computationally) to 
evaluate. The finite element concept may be applied to any of a variety 
of modeling problems; this investigation centered specifically on fine 


structure gravity modeling. 


Previously, earth-fixed spherical coordinate based finite element 
fields have been successfully used to replace globally valid spherical 
harmonic representations of the geopotential and its derivatives. But 
because the model used to simulate gravity disturbance data (Model 310) 
consists of a set of point masses placed in relation to the Geodetic 
Reference Surface of 1967, geodetic coordinates (H, 4, ¢ - see Figure 1) 


were used in the finite element modeling process. 


Finite element fields may be defined with geodetic coordinates as 
shown in Figure 2. The region to be modeled is defined by a set of 
maximum bounds (Hoax! Ne? Vase? and a set of minimum bounds (Hain! Amin’ 
Simin)’ The region is then divided into a set of smaller volumes, or 
finite element cells, where each cell has dimensions Heell’ dcell’ % sell" 
Thus the 300 km x 10° x 10° region above the reference ellipsoid, 

altitude 0 =- 300 km, 
longitude 70°E - 80°E, and 
latitude 25°S - 35°S, 


could be modeled by one hundred (100) cells each 300 km x 1° x 1° or 


four (4) cells each 300 km x 5° x 5° or one cell 300 km x 10° x 10°, etc. 
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earth-fixed unit vectors 


cos $ cos \ cos $¢ sin A sin $ i 
= - sin A cos 4 ie} j 
- sin $¢ cos \ -sindsindA cos > k 


Figure 1 Earth-Fixed Rectangular 
Coordinate Systems. 


instantaneous up, east, and north unit vectors. 


and Geodetic 


Figure 2 Geodetic F 


The approximation valid in any one cell is generated independentally 
of the approximation for any other cell. A set of gravity disturbance 

| "observations" are produced for each cell directly from Model 310 and 

then approximated for only that cell in a least squares fitting process. 

The approximating models thus generated are then stored on a high speed 

rotating mass storage device for later use in the computation of the 

approximate gravity disturbance in any finite element cell of the region 


modeled. 


ly OF Ak es mr Fe * : * a Re 


ERTL M0 Ai a TET knew 


4.0 FINITE ELEMENT FORMULATION 
Define the gravity disturbance, 6g: 
Zactual ~ Lreference * ég 
8g 


up 


ag: OF nat * Sactual ~ Zreference 


Oo een 


Now, replace the components of 6g with locally valid Chebyshev polynomial 
approximations, 6G, of the form, 


¥, (817, O14, Gy) 


(3) 
where, 
NORDER = the highest order Chebyshev polynomials, 


m=2- 3 


(ho H eg) /H 


= (x = deep) 


cell 


cell 


= = Wage) Pett 


TX) shifted Chebyshev polynomials, and 


Cups 5x’ Ceasti 5," and Cnorthi4k are constant ea determined via 
least square fits so that y (8G 2 y (6G 


up Zi ee 


- 6g ),*, 
i=l in 


east i 


m 
and ) (6G 
i=1 


a 2 ae ee 
— 69 orth) i are minimized for some local volume of space 
(i.e. the cell being modeled). 


Chebyshev polynomials, ti (x), of order n may be generated: 


SE I PT ay, 


See oe 


n t(%) 
fe) } 

1 Xx 

2 2n* - 1 


For n > 2, a recursion may be used: 
t (4) = 2x t) (x) - tk), ~“E< e <4 (5) 
Shifted Chebyshev polynomials, T(x) may be computed by substituting 


X= 2x -1 
so that, 
TX) 


t (X) = t (2X - 1), 6<z% <1 (6) 
n n — = 


To generate shifted Chebyshev polynomials a set of X's with values 
between 0 and 1 are needed. The difference between some point (H,A,$) 

in a cell and a reference point (HrefrAre¢: bref) (the “lowermost corner" 
of the cell) is divided by the cell's dimensions, Hoej1r Acell: %cell: 
for 


as in (4) to obtain a set of non-dimensional coordinates, Xie Kg 


2 BF 
each cell. The non-dimensional coordinates may then be used to generate 


the shifted Chebyshev polynomials. 


To determine the three sets of constant coefficients, Cups Coasts 


and Crorth’ for each cell, formulate and solve the linear least squares 
problems, 
[186,,,- S4ypI = [1ACy, ~ 94g) = min. 
[18e..ce - Seas! ! = [aC sce - SSeace! | = min- nae 
[18¢, rth io 89 orth! | is [lac orth rs 89. orth! | ——— 


where the coefficient matrix, A, has the form, 


7, (X3,) 1, (%2,) 7, (%,,) 7X3), (X3,)T, (RQ, 
By Mgq! Ty Sag? ag? eae) ae 


. e . ° . e 
e ° e ° . . 
e . . . . . 


T, (X30) 7, (XQ) TL (Xp) TS (Xa) (Xy Ty, (Xa) oe 


1 
1 
(8) 
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and the three coefficient vectors, C ,C 7, ana C , have the form 
| up east north 
(a; € ec 
YPo009 eastoog north 4 
Cc te. Cc 
UP O01 eastyo) north) 9, 
Cc Cc (e 
UPo10 eastyig northsi9 
(S Cc c 
} = up = east eS north 
| eS = 100 et 100 . Ce 100 (9) 
€ 5 Cc 
UPo92 eastyo> north, 95 
c c 
UP O11 eastyiy north), 
€ Ec 
UPr00 eastioo north, 49 
where 
N = NORDER, 
and the three sets of gravity disturbance observations, Sa OF nat" and 
S9,orth’ = 
&g (Xie oe) ég (X pee Gh) 
RE | up 1, 21 34 east 1, 21 35 
ég (X re. 7X ) ég (X Xx Pps ) 
| up 1, 25 35 east 1, 2, 3, 
BE i 6g = Sg a 
; up 6g 2 4X, 7X, ) east ég (x. 7X5. 7%, ) 
E up Ae. 23 33 east 73 23 33 
' a e . 
F j . ° 
i un 1 fie ie ) "Seane ae Mae ) 
f m m m m m m 
f 6 (X X x. ) 
H Snorth ae ae 
i (10) 
fr 6g (X 1X 7X ) 
| north 1, 25 35 
6g = 
north 89 orth (*1 X, X, ) 
. 2 3 o 
So north "1 1% ir ) j 
m m m 
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where {(X 1% Xe 
ZR i 
coordinates - usually a uniform observation grid in any one cell. The 


1 y, 1 =1, 2, ..., m} = a specified set of (Hy A515) 


gravity disturbance observations are determined directly from Model 310 


evaluations on the observation grid in each cell. 
4.1 Discussion 


Certain aspects of the process by which finite element fields are 
created may be capitalized on under special circumstances. When it is 
known beforehand that uniformly gridded gravity disturbance observations 
will be available, only one least squares matrix (A matrix) need be 
generated. This is true because the set of non-dimensional coordinates, 
COhy hy ey, 


5 i ¥ ; : ; : : 
the same as @very other cell's non-dimensional observation grid coordinates 


X,), i = 1, 2; ..., m} of one cell's observation grid will be 
when the observation grid pattern is set for all cells. If the same 
number of observations are made in each cell and the positions of the 
observations points are the same relative to the reference point (HY 
r 


ef’ 
, o ) in each cell, then all cells will have the same A matrix. 

ref ref 

The A matrix, after all, consists only of products of Chebyshev polynomials, 

which are functions only of the non-dimensional coordinates of the obser- 


vation points. 


This all means that only one A-matrix need be generated for the 
first cell of a finite element field; but the gravity observations must be 
generated for each cell. A reduction of the A-matrix to upper triangular 
form may be performed for the least Squares fitting process, after which 
each set of observations may be similarly reduced and then back-substituted 
to produce a set of coefficients for each cell. The A-matrix, again, 


needs to be reduced only once. 


It should be noted that this method of finite element field generation 
does not rely on Model 310 or any other particular gravity disturbance 
model. Model 310 could be replaced by any other process that is capable 
of producing the gravity disturbance observations at the observation 


grid points of each cell. 


Once a finite element field has been created, the Chebyshev polynomial 


gravity disturbance approximation 6G may be computed by evaluating (3). 
Lz 
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The subscripts ijk always follow the same pattern for a given NORDER. 


This means that a “current coefficient number", &, may be attached to 
each group of subscripts ijk and that the total number of coefficients, 
NC, necessary to model any one component of the gravity disturbance in 
one cell will be constant for a given NORDER. Table 1 shows these 


relationships for various NORDER's. 


Since we can precompute and save the subscripts, ijk, in the appro- 
priate pattern such that they are functions of only the current coefficient 


number, £, (3) can be written, 


(S 
6G up 
up Ne £ 
= 6G = 
a east Ceast ee tae is ae = 
2=1 kL 
ae cen c 
northp 
! 
| where, 
i NC = the total number of coefficients for one component of the 
gravity disturbance, as above 
ijk = precomputed and saved as functions of 2 
XX 0X, = functions of H,A,¢? and the constants Het’ det’ a 
H as in (4) above. The total number of summations in 


cell’ *ce11’ *ce12 
(11) is the same as in (3) but the amount of bookkeeping has been cut 


down tremendously by doing business in this fashion. Computational 
i savings of up to 50% were obtained by using (11) in lieu of (3) for 
evaluation purposes. Table 1 was also used in the generation of the A 
matrix, but the saving here is not as significant due to the previously 


discussed fact that the A matrix is only generated once. 


3 
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NORDER ijk 2 oa 
00 


i 001 2 
010 
100 4 4 
2 002 5 
OLX 6 
020 7 
101 8 
L160 9 
2°00 10 10 
3 00 3 ll 
Or 2 12 
N OON . 
0 1 N-1 . 
eo ce, ° N+1 n 
NOO NC ) ) n 
n=1\ m=1 


Table 1 Subscript Patterns for Various NORDER's. 
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5.0 DISCUSSION OF SOFTWARE 


A variety of FORTRAN programs and subroutines were written during 

the course of this research effort. The creation of finite element 
gravity disturbance fields may be accomplished by the execution of PROGRAM 
LOCALG followed by the execution of PROGRAM FINEG. Once a finite element 
field has been generated, the gravity disturbance approximation may be 
calculated at any point in the finite element field by SUBROUTINE FINITE. 
To calculate the gravity disturbance from the mass points of Model 310, 
SUBROUTINE PTMASS was written. SUBROUTINE PTMASS was used both to produce 
observations for the least squares fitting process and to produce obser- 

{ vations to test against for error analyses of finite element fields. 
Because SUBROUTINE PTMASS required the precomputed products of each 
point mass and the gravitational constant and the earth-fixed Cartesian 
coordinates of each point mass of Model 310, another small program, 


PROGRAM MASPOS, was written to calculate and store these quantities. 


PROGRAM LOCALG requires input from the user which defines the finite 
| element field to be generated. It writes the coefficient matrix (A 
matrix) and a set of gravity disturbance observations for each cell of 
the finite element field onto a magnetic tape (TAPE 3 - this may be any 
sequential mass storage device) for later use by PROGRAM FINEG. PROGRAM 


LOCALG is the first part of the local gravity disturbance modeling process. 


Figure 3 shows the logic flow for PROGRAM LOCALG. 


The second part of the local gravity disturbance modeling process - 
the reduction of the A matrix and the calculation of the fine structure 


gravity model coefficients (9) - is accomplished by PROGRAM FINEG. 


‘ 

4 

i 

; 

i 

j PROGRAM FINEG reads the A matrix and observation from the magnetic tape 

(really, sequential file) created by PROGRAM LOCALG. PROGRAM FINEG 

i reduces the A matrix and calculates each cell's model coefficients in a 

i least squares reduction. The model coefficients are written onto a 
random access file (preferably a high speed disk file but any high speed 
mass storage random access file could be used). The logic flow for 
PROGRAM FINEG is shown in Figure 4. 
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START 


initialize constants 
! Ore 
read point mass values and precomputed positions qv 


input Chebyshev polynomial order, observation grid pattern QuVWwY 7 
finite element cell size, and finite element field bounds NORDER, MOBSU, MOBSE, MOBSN. 
HCELL, ALCELL, APCELL 


compute indexes, indicators, and increments for foe ALMAX? ae 


finite element field 


write critical finite element field information onto the VW AAA> ( Tape 3 
sequential file 


for one cell of the finite element field loop through 
the cell, one layer of the observation grid (H=constant) 
per loop (total of MOBSU loops) 


generate Chebyshev polynomials on the observation grid 
and store in the A matrix 


write the A matrix onto the sequential file (one wwe ( Tape 3 
layer will be written each loop) 


next layer in cell 
loop through finite element field, one cell per loop 


generate gravity anomaly “observations,” sg , 
(CALL PTMASS) i: 


write observations onto the sequential file wwe (0) Tape 3 


next cell in finite element field 


59s 59, 


STOP 


Figure 3 PROGRAM LOCALG Logic Flow. 
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START 


j Tape 3 
read critical finite element field information from «m(0) 


the sequential file 


write critical finite element field information onto the VVVVp> 
first record of the random access file 


loop through one cell of the finite element field, 


one layer per loop 


the A matrix 


next layer 


reduce the least squares matrix to upper triangular 
form (CALL ALSQ) 


loop through the finite element field, one cell 
per loop 


from the sequential file 


compute polynomial coefficients c of the finite 


element model by minimizing || Ac-sq obs. || 


write coefficients of the model for this cell onto YANN 
the random access file 


next cell 


STOP 


Figure 4 PROGRAM FINEG Logic Flow. 


read one layer of the least squares matrix fron </VWVVVS (D)rre 3 
the sequential file into the appropriate locations in 


read gravity anomaly observations for this cell qe“ (un 3 
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SUBROUTINE FINITE evaluates (11) above at some point (H,A,$). 
Using the sets of finite element coefficients previously generated by 
PROGRAM FINEG, SUBROUTINE FINITE returns an evaluation of the gravity 
disturbance approximation (11) at the specified point. The logic flow 


for SUBROUTINE FINITE is shown in Figure 5. 


Gravity disturbance observations were produced for this study by 
SUBROUTINE PTMASS, which uses the 1080 point masses of Model 310 and 


(12) below to compute the gravity anomaly, 65g, at some point P(x,y,z): 


6 
Sy 1080 km, 
6 = 6 = _—_ * 
g ay ) 5 Ay (12) 
ég i=l da; 
z 
where 
k = 6.67 x 10714 m3/(g+sec*) = the gravitational constant 
m =1x 10!9, 0, or 1 x 10!9g 
da x- x 
x. i 
i 
(13) 
-i i 
d zZ-2Z 
Zz. i 
i 
d,=/vd 24d 2+4d._.2 
i x. y. Zz. 
i i i 


To minimize the execution time of SUBROUTINE PTMASS, PROGRAM MASPOS 
was written to compute and save the products of the gravitational constant 
and each of the 1080 point masses (14) and the earth-fixed x,y,z coordi- 


nates of each of the 1080 point masses (15): 


(km),, i = 1, 2, ..., 1080 (14) 
Xin Yur Ze i= 1, 2; «22; LOSO (15) 
18 
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CALL FINITE (specify H,A,@ coordinates of some point) 


Disk 


ical finite el co 
read critica nite element 
field information from the www 
first record of the Tape 1 
random access file 


point in 
finite element 
field? 


cal] PTMASS 


RETURN 


compute indexes 


determine which cell of the finite element field the 
point lies in (calculate IPOINT, the cell pointer) 


isk 


read appropriate mow 
coefficients Tape 1 


coefficients 
of appropriate cell 
in core? 


compute finite element cell boundaries; non-dimensional ize 
the calling coordinates; generate Chebyshev polynomials 


Compute 5G 
(matrix product of finite element model coefficients and 
Chebyshev polynomials) 


RETURN 


Figure 5 SUBROUTINE FINITE Logic Flow. 
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6.0 TRADEOFF STUDIES 


Several tradeoff studies were conducted on the UVA CDC Cyber 172. 
The central processor execution times of SUBROUTINE FINITE were compared 
to SUBROUTINE PTMASS (Model 310). Although it is really a price that is 
paid only once and a priori in the laboratory, the central processor 
execution times necessary to create finite element fields were examined. 
Consideration was given to the tradeoff: total number of finite element 
coefficients versus NORDER versus cell size, for finite element models 
of a given region. Lastly, error analyses for a variety of NORDER's and 
cell sizes were conducted; specifically, the maximum absolute error versus 


NORDER versus cell size tradeoff was examined. 


Figure 6 shows the execution time test results for SUBROUTINE FINITE 
and SUBROUTINE PTMASS. SUBROUTINE PTMASS has a nearly constant run time. 
The run time for SUBROUTINE FINITE, however, depends on NORDER, the 
Chebyshev polynomial order, and on whether or not a random access call 
is needed. A random access call is needed anytime the set of coefficients 
currently in core is not the set for the cell in which the calling 
coordinates lie; that is, the previous calling coordinates were ina 


different finite element cell than the current calling coordinates. 


Observe that the random access time is approximately constant (0.7 
ms) regardless of NORDER. This is true because most of the random access 
time is spent in finding the address of the first coefficient on the 
rotating mass storage device. The read time for one cell's set of 
coefficients (even NORDER = 6 has only 252 total coefficients for one 
cell) will be very small compared to the first word address find time. 
The execution time ratios show that even for the worst cases SUBROUTINE 
FINITE will be at least 20-50 times faster than an explicit point mass 


disturbance acceleration model (e.g., SUBROUTINE PTMASS). 


A certain price must be paid to generate finite element fields. 
Execution times involved in running PROGRAM LOCALG and PROGRAM FINEG are 
shown in Figure 7 for the UVA CDC Cyber 172. The total price for running 
both programs is also shown. The generation and reduction of the A 


matrix is a one-time price for any one finite element field. The price 
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SUBROUTINE PTMASS 
SUBROUTINE FINITE 


NORDER 


Nn nO B® WwW 


(SUBROUTINE FINITE Execution Time) 
(SUBROUTINE PTMASS Execution Time) 


NORDER 


3 
4 
5 
6 


Figure 6 


: 123.2 ms 
Execution Times (ms) 
No Random Access Random Access 
eT 2.4 
2.6 3.3 
3.8 " 
Se Gk 


(computational speed advantage) 


No Random Access Random Access 
1/72'.5 L/Si.s 
1/47.4 L/37.3 
1/32.4 1/27.4 
1/22.8 1/20.2 


SUBROUTINE FINITE vs. SUBROUTINE PTMASS 
Central Processor Execution Time Comparison. 
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NORDER 


3 4 5 6 
PROGRAM LOCALG 
time to generate A matrix (sec.) we «3 oe Tse 
time to generate observations 
per cell (sec.) 7.9 5.4 <26.6 42.3 
PROGRAM FINEG 
time to reduce A matrix (sec.) 4 20 .7 30.4 
time to calculate coefficients 
per cell (sec.) a2 «J So 4.3 ' 
TOTAL PRICES 
ij 
| A matrix and reduction (pay once) 1 
(sec.) .6 203 9.4 31.8 
observations and coefficients 
per cell (sec.) Sol 16.1 28.5 46.6 


Given: H = constant = 300 km 14 
cell j 


Figure 7 Central Processor Times to Create 
Finite Element Fields. 
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per cell must be multiplied by the total number of cells in the finite 
element field and added to the price to generate and reduce the A matrix. 
Note that over 90% of the price paid per cell is just the generation of the 
| gravity observations via SUBROUTINE PTMASS. 
f The total numbers of coefficients needed for various finite element 
: field models of the 35°x40° region covered by Model 310 are shown in Figure 
| 8. Figure 8 assumes that Heell will be constant at 300 km and that 300 km 
will be the maximum altitude to be modeled. For a given cell size, it 
f shows how many finite element cells would be needed to model the region 
and how many total coefficients would be needed for various NORDER's to 
model the region. Based upon error analyses conducted on a small volume 
near the center of the 35°x40° region, expected maximum error bounds were 
determined for certain cell sizes and NORDER's. The upper sol 
dashed underline indicates an expected maximum absolute error o# 
mgals and the lower solid and dashed underline indicates an expected 
maximum absolute error of 3.0 mgals. Solid underlines in both cases 
indicate that the expected maximum absolute error was verified for small 


test volumes while dashes indicate the apparent trend. 


In order to study the maximum absolute error as functions of NORDER 
and cell size, it was decided to consider the maximum absolute error as 
| a function of cell size only for fixed NORDER's (see Figure 9) and to 
| consider absolute error as a function of NORDER for fixed cell sizes 
(see Figure 10). For all cases, cell altitude (Hoeiy) was fixed at 300 


km. 


Figure 9 shows that the maximum absolute error varies approximately 
as the square of the cell dimension in A or $. For constant cell altitude 
and for small oell and $¢, the volume will vary approximately as the 
square of i or > if \ se] = doell' Heuce the maximum absolute 


cell cell’ 


| 

$ 

| 

oT 

| 

i 

| error is really proportional to cell volume. Figure 10 shows that the 
maximum absolute error varies approximately inversely as the square of 

NORDER. As the number of coefficients per cell varies approximately as 
the square of NORDER, the maximum absolute error may be considered to 


vary nearly inversely with the number of coefficients per cell. 


23 


x 


ee 
Bes fiteste de Viea.2 


TES TSS = Re ek xa ; OS at ae 


NORDER 
Cell Size 
Qoei1 x del? No. of Cells 3 4 S 6 
\ . 
\ \ 
\ \ 
1.0°x1.0° 35x40 = 1400 \, 84,000 \\ 147,000 235,200 352,800 
aa aaa ‘N a 
Fa25%i [25° 28x32 = 896 53,760 ‘\_ 94,080, 150,528 225,792 
pee Race es Tea aegis So 
~~ Ns 
P.S°xy. 5° 24x27 = 648 38,880 68,040 ».108,864 “y, 263.296 
%, 
Le 75exE. 75° 20x23 = 460 27,600 48,300 77,280 ot ETS, 920 
% 
\ 
20S x20 18x20 = 36 r 7,8 , 480 ..90,720, 
2,022.0 8x20 (0) 21,600 37,800 60 , + 
2.25 °X2 25° 16x18 = 288 17,286 30,240 48,384 72,576 
25°? 5° 14x16 = 224 13,440 23,520 37,032 56,448 
2752x207 5> 13x15 = 195 11,700 20,475 32,760 49,140 
302X307 12x14 = 168 10,080 17,640 28,224 42,336 
Given: H = constant = 300 km 
cell 


* ~ Finite element fields above this line are expected to have 
|exxor| wax < 1.5 mgals. 


‘ + ~ Finite element fields above this line are expected to have 
error | pax < 3.0 mgals. 


Figure 8 Total Number of Coefficients vs. NORDER vs. Cell Size 
For a 300 km x 35° x 40° Finite Element Field. 
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Figure 9 Maximum Absolute Error versus Cell Size 
for Fixed NORDER's. 
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a. Cell Size = 1°x1°x300 km b. Cell Size = 1.5°*1.5°*300 km 


(mgals) 


c. Cell Size = 2°x2°x300 km d. Cell Size = 3°x3°x300 km 
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NORDER 


Figure 10 Maximum Absolute Error versus NORDER 
for Fixed Cell Sizes. 
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7.0 EXAMPLE FINITE ELEMENT GRAVITY MODELS 


Based upon the tradeoff studies conducted as described above, two 
example finite element gravity fields were generated. The 300 km x 10° x 
10° region with altitude 0 to 300 km, longitude 70°E to 80°E, and latitude 
25°S to 35°S was modeled. A finite element field with NORDER = 3 and 
with one hundred 300 km x 1° x 1° cells was generated as was a finite 
element field with NORDER = 5 and the same size and number of cells. An 
error analysis was conducted on both fields. The finite element model 
was compared on a point by point basis with Model 310 at sample points 
“away from the fit points" for altitudes of 1m, 150 km, and 299 km. The 
error mean, the RMS error, and the maximum absolute error were determined 
for the three sets of sample points for both finite element fields and 


are shown in Figure 11 and Figure 12. 


6g 6g 


up east north 
(mgals) (mgals) (mgals) 
H=1mMm 
error -.005 -.001 004 
RMS error % -963 -645 Pay | 
[error] wax 2.717 2.074 2.236 
H_= 150 km 
error -.023 -003 -009 
RMS error «252 224 270 
| Jerror| wax 1.149 .701 1.039 
H = 299 km 
error -018 -.006 -.007 
| RMS error 424 .277 . 322 | 
|extor | max 1.130 ae 1.002 | 


tee me eT 


Figure 11 NORDER = 3 Example Finite Element Gravity 
Disturbance Field. 
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(mgals) (mgals) (mgals) 
B= lm 
error -.004 -001 -002 
RMS error 224 ~ 150 -167 | 
|erxor| wax 918 637 581 | 
| H = 150 km | 
error .007 -.002 -.004 | 
RMS error .O71 .044 -055 | 
| |ezrer | ax 216 als? a2l5 | 
| 
H = 299 km 
error -.002 001 001 
| RMS error 101 067 JO77 
[error | ax .329 254 249 


Disturbance Field. ) 
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| 
| 
Figure 12 NORDER = 5 Example Finite Element Gravity 
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8.0 CONCLUSIONS 


This effort has demonstrated that: 

« Finite element fields are suitable for modeling the gravity 
disturbance. 

* Finite element fields are computationally attractive for modeling 


the gravity disturbance. 
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APPENDIX - SOFTWARE 


A.OQ FORTRAN IV programs have been developed which create finite element 
gravity disturbance fields and which evaluate the gravity disturbance 
within these fields via the methods discussed in this report. Descriptions 
of these routines and their attendant subroutines are provided here, 


together with listings of all routines and example executions. 


Three major routines were developed - PROGRAM LOCALG, PROGRAM FINEG, 
and SUBROUTINE FINITE - to create finite element gravity disturbance 
fields and to evaluate the gravity disturbance within these fields. Two 
additional routines were developed to provide simulated gravity disturbance 
data - PROGRAM MASPOS, which, for the 1080 mass points of Model 310, 
computes and stores the products of each mass point and the gravitational 
constant together with the geocentric rectangular coordinates of each 
Mass point, and SUBROUTINE PTMASS, which evaluates the gravity disturbance 
as Simulated by Model 310. 


A.1.1 PROGRAM LOCALG generates the A matrix (8) for a "typical" finite 
element cell and a set of gravity disturbance observations (10) for each 
cell of a finite element gravity disturbance field. The gravity distur- 
bance observations are evaluations of SUBROUTINE PTMASS on the uniform 
observation grid in each cell. The user defines a finite element gravity 
disturbance field by providing the following items as inputs to PROGRAM 
LOCALG: 

1. The integral order of the locally valid Chebyshev polynomial 
approximating functions, NORDER, (0 < NORDER < 6). 

2. The observation grid pattern (same for all cells) to include 
the number of observations in the up, eastern and northern directions, 
MOBSU, MOBSE, MOBSN. The total number of observations (MOBSU x MOBSE 
x MOBSN) should be at least three times the number of coefficients, 
NC (see Table 1), in the approximating model (3). Furthermore, 
the number of observations in each direction should be equal (MOBSU = 
MOBSE = MOBSN). A good rule of thumb is that for all NORDER's 


between 3 and 6 inclusive, the number of observations in each direction 


should be one more than NORDER (MOBSU MOBSE = MOBSN NORDER + 1). 


2) 


The maximum number of observations allowed by the program is 343 
| (a 7X7x7 observation grid). Non-uniform sample grids (e.g. 5x7x7, 
6x5x4, etc.) produce unpredictable results - often an unacceptable 
fit. 
3. The size of each finite element cell in H,A,¢ must be specified, 
HCELL, ALCELL, APCELL. HCELL is cell size in H (meters), (HCELL < 
{[HMAX - HMIN]). ALCELL is cell size in A (degrees), (ALCELL < {ALMAX 
- ALMIN]). APCELL is cell size in $ (degrees), (APCELL < [APMAX - 
APMIN]). 
4. The minimum and maximum H,A,?> boundaries of the finite element 
field, HMIN, ALMIN, APMIN, HMAX, ALMAX, and APMAX. HMIN and HMAX 
are heights (H) above the reference ellipsoid in meters, (HMIN > 0 
and HMAX > HMIN). ALMIN and ALMAX are geodetic longitudes (A) in 
degrees, (-180° < ALMIN < ALMAX < 180°). APMIN and APMAX are 
geodetic latitudes ($) in degrees, (-90°< APMIN < APMAX < 90°). 


| These thirteen items are input on four data cards. NORDER, MOBSU, 
MOBSE, and MOBSN are input in (4110) format on the first card. The 
remaining three cards must have a (3E20.14) format, with the 
finite element field delimiters, HCELL, ALCELL, and APCELL on the second 
card, HMIN, ALMIN, and APMIN on the third card, and HMAX, ALMAX, and 
APMAX on the last card. 


j Because gravity disturhance observations are generated within 
PROGRAM LOCALG by SUBROUTINE PTMASS, PROGRAM LOCALG also requires as 

H input the precomputed products of the gravitational constant and each of 
the mass points of Model 310 and the precomputed rectangular coordinates 
of each mass point. These quantities should have previously been stored 
sequentially via an unformatted write on TAPE 2. PROGRAM LOCALG accesses 


these quantities with an unformatted read of TAPE 2. 


PROGRAM LOCALG produces as output an unformatted sequential file, 
TAPE 3, which contains the A matrix for a typical finite element cell 
and a set of gravity disturbance observations for each cell of the user 
specified finite element field. This file is properly formatted for 
later use by PROGRAM FINEG. PROGRAM LOCALG requires two supporting 
subroutines - SUBROUTINE CHEBY and SUBROUTINE PTMASS. 
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A.1.2 SUBROUTINE CHEBY (X, N, TA) returns Chebyshev polynomials in array 
TA, evaluated at a specified value, X, through the specified order, N. 
TA(1) = Ty (X), TA(2) = T(X). TA(3) = T, (xX), etc. TA is a seven element 


array; hence the specified order is bounded, 0 < N < 6 


A.1.3 SUBROUTINE PTMASS produces gravity disturbance “observations", 


and ég from (10) to be used by PROGRAM LOCALG and 


a SJeast’ north’ 
PROGRAM FINEG for fitting purposes and to be used to test as data against 
SUBROUTINE FINITE. SUBROUTINE PTMASS produces the gravity disturbance 
data by using the 1080 point masses of Model 310. SUBROUTINE PTMASS 
requires that COMMON /MASPOS/ be filled with the precomputed products 

of the gravitational constant and each point mass and the geocentric 
rectangular coordinates of each point mass. SUBROUTINE PTMASS interfaces 
with the calling program via COMMON /XYZ/ through which the rectangular 
coordinates, x, y, z, of the calling point are received and the components 
of the gravity disturbance Sg. aa and ég., (DELGX, DELGY, DELGZ) are 
transmitted. The rectangular components of the gravity disturbance must 
and ég 


be transformed to geodetic components 6g h’ outside 


6 
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SUBROUTINE PTMASS (i.e., in the calling routine). 


A.2.1 PROGRAM FINEG generates sets of locally valid coefficients epi 
c and c. (9) for each cell of a finite element field which was 


east’ orth 
previously specified by inputs to PROGRAM LOCALG. The only input required 
by PROGRAM FINEG is the unformatted sequential file, TAPE 3, previously 
produced by PROGRAM LOCALG. PROGRAM FINEG produces as output a random 
access file, TAPE 1, which contains, in effect, a finite element gravity 
disturbance field. The first record of the random access file contains 
certain important information about the finite element field (NORDER, 
boundaries, etc.). In addition to the first record, the random access 
file consists of one record of coefficients (9) for each cell of the 
finite element field. Thus the number of records in the random access 
file will be equal to the number of cells in the finite element field 
plus 1. PROGRAM FINEG requires one supporting subroutine - SUBROUTINE 


ALSQ. 


A.2.2 SUBROUTINE ALSQ (A,Y,B,R2,NN,MM,NA) solves the linear least squares 
problem, | |aB = y|| = minimum. SUBROUTINE ALSQ requires as input the 
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coefficient matrix, A, the observation vector which is to be fit, Y, the 


number of rows used in the A matrix, NN, the number of columns used in 
the A matrix, MM, and the first dimension of the A matrix, NA. NA must 
be equal to one plus the maximum number of rows in the A matrix (NA = 
NN ax + 1). SUBROUTINE ALSQ returns B, the coefficients of the fit and 
R2, the sum of the squares of the residuals. 

A.3.1 PROGRAM FINTES is a simple error analysis routine. It compares 
the gravity disturbance as approximated by SUBROUTINE FINITE with the | 
gravity disturbance as evaluated by SUBROUTINE PTMASS (Model 310). 

Residuals are printed for each point of a user specified grid, along 

with the mean of the residuals, the root mean square of the residuals, 

and the maximum absolute value of the residuals, for each component of 
the gravity disturbance. The following items are required as input to | 
specify the sample grid: 

l. The number of samples in each direction, ISTEPH, ISTEPL, and 
ISTEPP, (ISTEPH > 0, ISTEPL > 0, and ISTEPP > 0). The total number 
of sample points on the grid will be the product of these three 
numbers. 

2. The minimum sample grid coordinate values ("the lower corner"), 
HMIN (meters), ALMIN (degrees), and APMIN (degrees) and the maximum 
sample grid coordinate values ("the upper corner"), HMAX (meters), 
ALMAX (degrees), and APMAX (degrees), (0 < HMIN < HMAX, ~180° < 
ALMIN < ALMAX < 180°, -90° < APMIN < APMAX < 90°). 

These nine items are input on three data cards. ISTEPH, ISTEPL, and 
ISTEPP are input on the first data card in a (315) format. HMIN, ALMIN, 
and APMIN are input on the Second card in (F10.0, 2F10.2) format and 
HMAX, ALMAX, and APMAX are input on the third card in the same format. 
PROGRAM FINTES requires three supporting subroutines - SUBROUTINE FINITE, 
SUBROUTINE CHEBY (see discussion of this routine above), and SUBROUTINE 
PTMASS (also discussed above). 


Because PROGRAM FINTES uses SUBROUTINE PTMASS, it also requires as 
input the unformatted sequential file, TAPE 2, containing the precomputed 
products of the gravitational constant and each point mass of Model 310 


and the precomputed rectangular coordinates of each point mass. 
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4.3.2 SUBROUTINE FINITE directly replaces Model 310 in applications 
involving the gravity disturbance 6g. Given the geodetic coordinates, H, 
ALAM, and APHI (H,A,$), of a specified point, SUBROUTINE FINITE returns 


the geodetic components, GU, GE, and GN (6G, , 5G ) of the 


Pp east’ "6 oeths 
gravity disturbance approximation, 5G (3). SUBROUTINE FINITE determines 
which set of finite element coefficients to use at the specified point, 
reads them into core, and computes 6G. SUBROUTINE FINITE requires 
SUBROUTINE CHERY (described previously) to produce Chebyshev polynomials 
at the specified point. Whenever the specified point does not lie within 
the bounds of the finite element field currently accessed, SUBROUTINE 
FINITE calls SUBROUTINE PTMASS, writes an error meSsage, and sets the 


error flag ISITIN to zero. 


Prior to the first call to SUBROUTINE FINITE, the random access file, 
Tape 1, must be made available to the calling program, and two flags, 
IFLAG and ISITIN, must be set to zero and included in a COMMON in both 
the calling program and SUBROUTINE FINITE. This could be COMMON /HLP/ 
if necessary (COMMON /HLP/ is the vehicle by which SUBROUTINE FINITE 
receives the calling coordinates and returns the components of 6G), 


although COMMON /IMARK/ was used. 


The requirements of CDC FORTRAN dictated that the CDC utility 
subroutines, OPENMS, WRITMS, READMS, and CLOSMS, be used to access the 
random access file. Furthermore, these CDC utility routines require that 
the user establish an index array (IMARK, in this case) for the random 
access file. If an index array is required by a particular version of 
FORTRAN, the index array, IMARK, should be included in a COMMON in both 
the calling program and SUBROUTINE FINITE. 


A.4 PROGRAM MASPOS creates an unformatted sequential file, TAPE 2, which 
contains the precomputed products of the gravitational constant and 

each point mass of Model 310 and the precomputed geocentric rectangular 
coordinates of each point mass. PROGRAM MASPOS requires as input coded 
information about Model 310 in the form of 30 cards. Each card represents 
36 mass points in a grid row of equal latitude where the mass points are 
50 longitude minutes apart. The mass points should be listed on each 


card by increasing longitude and the cards should be arranged in order 
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of decreasing latitude. PROGRAM MASPOS reads only the coded multiplication 


factors for each of the 36 mass points of the grid row, where the coded 


factors indicate the mass of the respective mass points, 1 = -1.x10!9 grams, 


2 - 0. grams, and 3 = +1.x1019 grams. The coded factors are located in 


columns 25-60 of each data card. All other data card columns are ignored 


(for information only, columns 1-9 indicate south latituce in arc minutes 


of the grid row and columns 10-18 indicate east longitude in arc minutes 


of the west-most mass point in the row). 
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